//============================================================================//
//                                                                            //
//                              Program/Module                                //
//                                   from                                     //
//                 Selfverifying Solvers for Dense Systems of                 //
//                     Linear Equations Realized in C-XSC                     //
//                                                                            //
//                     Carlos Holbig and Walter Kraemer, 2003.                //
//                                                                            //
//       Program developed in C-XSC by Bernardo Frederes Kramer Alcalde,      //
//           Paulo Sergio Morandi Junior and Carlos Amaral Holbig             //
//                                                                            //
//============================================================================//
//----------------------------------------------------------------------------//
// File: lss_aprx (implementation)                                            //
// Purpose: Compute an approximate inverse of the square matrix A             //
// Global functions:                                                          //
//   MINV(): Computes the approximate inverse of a square matrix A using the  //
//           Gauss-Jordan algorithm.                                          //
//----------------------------------------------------------------------------//

#include "lss_aprx.hpp"

using namespace cxsc;
using namespace std;

//---------------------------------------------------------------------------
// minv computes approximate inverse of the matrix A by use of the
// Gauss-Jordan algorithm.
// A   : input matrix and outpout of approximate inverse
// errorCode : Error indicator = 0 everything ok, = 1 matrix singular
//---------------------------------------------------------------------------
int minv(rmatrix& A) {
	int       i, j, k, au1, au2, ao1, ao2, errorCode;
	bool      ok;
	real      h;
	rvector   aux( Lb(A,ROW),Ub(A,ROW) );
	intvector v( Lb(A,ROW),Ub(A,ROW) );
	int       l; 

	const real eps = 1e-15;

	au1 = Lb(A,ROW);
	ao1 = Ub(A,ROW);
	
	au2 = Lb(A,COL);
	ao2 = Ub(A,COL);

	ok = ( (ao1 - au1) == (ao2 - au2) )? true : false; // square matrix?

	// Vector that will control what's happening with the column
	for(i = au1; i <= ao1; i++){
		v[i] = i;
	}

	i = au1-1;  // Row index
	j = au2-1;  // Column index

	while ((j < ao2) && ok){
		i++; 
		j++;  // i and j run simultaniously

		// pivot search
		if (j < ao2){
			l = i;
			for(k = i+1; k <= ao1; k++)
				if(abs(A[k][j]) > abs(A[l][j]))
					l = k;

			// row interchange
			if(i != l){
				aux = A[l]; 
				A[l] = A[i]; 
				A[i] = aux;
				k = v[l]; 
				v[l] = v[i];
				v[i] = k;
			}
		}

		// transformation
		if(A[i][j] == 0){
			ok = false;
		} else {
			h = 1.0 / A[i][j];
			Col(A,j) = h * rvector( Col(A,j) );
			A[i][j] = h;
			
			for(k = au2; k <= ao2; k++){
				if(k != j){
					for(l = au1; l <= ao1; l++){
					    if(l != i){
							h = A[l][k] - A[l][j] * A[i][k];
							if(h == 0.0) {
								A[l][k] *= eps;
							} else {
								A[l][k] = h;
							}
						}
					}
					A[i][k] *= -A[i][j];  // A[i][k] = -A[i][k] * A[i][j]
				}
			}
		}
	}

	// column interchange
	if(ok){
		for(k = au1; k <= ao1; k++){
			l = au2 + k - au1;
		    while(v[k] != k){
				i = v[k]; 
				v[k] = v[i]; 
				v[i] = i;
				j = au2 + i - au1;
				aux = rvector( Col(A,j) );
				Col(A,j) = Col(A,l);
				Col(A,l) = aux;
			}
		}
	}

	errorCode = !ok;
	
	return errorCode;
} // end MINV
